Detecting dispersion problems in GLMMs

Replicability crisis in Ecology too?

High false positive rates (type I error) in Ecological studies
A non-intentional cause is failing in checking model’s assumptions
Researchers don’t check models -> higher chances of false conclusions
Few tools for model diagnostics of GLMs, GLMMs.

Can you trust your model?

Dispersion problems in count data

Example count data:
- Species richness
- Abundance of individuals
Modeling with Poisson GLMs/GLMMs

UNDER or OVERDISPERSION:

When data has more or less variability than expected by the distribution used for modeling.

I’m pretty sure that many in the audience had modeled count data, for example, species richness or individual abundances. Here I’m bringing a very simple example of a field study evaluating if the abundance of a plant species changes with a environmental gradient, for example soil moisture. Those datasets are usually modeled using GLMs with a Poisson distribution. However, a problem that is very common in ecology is the phenomenon of over or underdispersion, that is, when the dataset has more or less variability than expected by the distribution you are using to model it. It happens because the Poisson distributions, for example, has a fixed mean-variance relationship: the mean of the distribution is also the variance. That is, if you have lower mean, you have lower variance, and vice-versa but it is not always the case in ecological data This problem is more common that we may expect, simply because we don’t check it properly.

GOALS

Aware ecologists of dispersion problems with count data

Identify and describe the 3 main causes by using model diagnostics tools with the DHARMa R package

Show some modeling solutions for these causes, with the glmmTMB R package

3 causes of dispersion problems

“Real” overdispersion:

Abundances vary more that expected by the model, in general.

Heteroscedasticity:

Zero-inflation:

Let’s start with the 3 main causes of having more or less variability in your data than expected by your model. Here I will focus on the overdipsersion problem, for simplicity, and use the example I showed before with the Poisson distribution. So we have the species abundances along an enviromental gradient. In green, we have the data behaving like a poisson distribution, so the species abundances increase with the gradient but there is no exta variability, then the fitted mean by the GLM model is the green line. In red, I have a the datasets with larger variation than expected by the Poisson model. We see that the red points are more spread thant the green one. This is what I’m calling here “real” overdispersion in contrast with the two other “apparent”overdispersion. Some would also call it Extra-Poisson variance and we could, for example, model,it with a Negative Binomial GLM.

3 causes of dispersion problems

“Real” overdispersion:

Abundances vary more that expected by the model, in general.

Heteroscedasticity:

Abundances variation increases with the environmental gradient.

Zero-inflation:

3 causes of dispersion problems

“Real” overdispersion:

Abundances vary more that expected by the model, in general.

Heteroscedasticity:

Abundances variation increases with the environmental gradient.

Zero-inflation:

More zero abundances than expected by the model.

Consequences of dispersion problems

OVERDISPERSION

Too small standard error / confidence intervals
Larger chance of false positive results

Wrong estimates, especialy if ignoring other processes (e.g. zero-inflation causes) in your data-generating process.

Missing the opportunity to learn more from your data. Ecological meanings for modeling unexpected variability.

If you ignore these problems in your model, the most clear consequence is that you have estimated too small standard errors for your estimates, this means that you think you have lower uncertainty, and smaller confidence intervals, but in reality your data says the opposite. This increases the chance of the type I error, that would tell you you have an significant effect, but this doesn’t exist. maybe that’s why people prefer to ignore these problems when fitting model to data ? Sometimes, depending on the complexity of your model and data, you may also have wrong estimates, that is, you estimate the slope of the realtionship of the abundance with your environmental variable, but this is wrong, and you will predict very poorly from your data. Especially if you have zero-inflation processes behind your data.

And what I really think is that, you are also ignoring how to learn more from your model, for example, if you find out that there is a variable that is increasing the variability in abundances, wouldn’t it be a cool ecological result?

Detecting dispersion problems with DHARMa

Residual diagnostics with `DHARMa`

Scaled quantile residuals -> Simulating from the model
Residuals between 0 and 1 for ANY model complexity or distribution
Interpreted the SAME way:

If your model is correctly specified, i.e. your have the “data-generating process”, scaled quantile residuals will present a uniform “flat” distribution between 0 and 1.

The DHARMa approach is to generate a specific type of residuals, the scale quantile residuals, by simulating many new datasets from your model. The advantange of these residuals is that it is standardized between 0 and 1 for any distribution and model complexity, and they should be interpreted the same way. That is, if your model is correctly specified, that means you have the data generating process, the residuals should present an uniform distribution, a flat distribution between zero and 1. These are the two main diagnostics plot from DHARMa, one quantile-quantile plot that you expect that the residuals (black dots) should align with the diagnoal red line, and the other the residuals agains the fitted values, which we would expect a uniform distribution of them along the fitted values, reprsented by these 3 black horizontal lines. Now, let’s check DHARMa residuals can detect dispersion problems in the data…

Detecting “real” overdispersion

Wrong model

m <- glmmTMB(observedResponse ~ Environment1 + (1|group),
        family = poisson(), data = overData)
res <- simulateResiduals(m)
plot(res)
testDispersion(res)

Dispersion = 5.19, p-value = 0.

Detecting “real” overdispersion

Wrong model

m <- glmmTMB(observedResponse ~ Environment1 + (1|group),
        family = poisson(), data = overData)
res <- simulateResiduals(m)
plot(res)
testDispersion(res)

Dispersion = 5.19, p-value = 0.

Solution

m <- glmmTMB(observedResponse ~ Environment1 + (1|group),
        family = nbinom2(), data = overData)
res <- simulateResiduals(m)
plot(res)
testDispersion(res)

Dispersion = 1.19, p-value = 0.224.

Detecting heteroscedasticity

Wrong model

m <- glmmTMB(observedResponse ~ Environment1 + (1|group),
        family = poisson(), data = overData)
res <- simulateResiduals(m)
plotResiduals(res, form = data$Environment1,
              absoluteDeviation = T)
testDispersion(res)

Dispersion = 1.9, p-value = 0.

Detecting heteroscedasticity

Wrong model

m <- glmmTMB(observedResponse ~ Environment1 + (1|group),
        family = poisson(), data = overData)
res <- simulateResiduals(m)
plotResiduals(res, form = data$Environment1,
              absoluteDeviation = T)
testDispersion(res)

Dispersion = 1.9, p-value = 0.

Solution

m <- glmmTMB(observedResponse ~ Environment1 + (1|group),
             dispformula = ~ Environment1, # dispersion formula
        family = nbinom2(), data = data) # but needs negative binomial
res <- simulateResiduals(m)
plotResiduals(res, form = data$Environment1,
              absoluteDeviation = T)
testDispersion(res)

Dispersion = 1.11, p-value = 0.44.

Detecting zero-inflation

Wrong model

m <- glmmTMB(observedResponse ~ Environment1 + (1|group),
        family = poisson(), data = overData)
res <- simulateResiduals(m)
plot(res)
testZeroInflation(res)

Zero-inflation = 5.16, p-value = 0.

Detecting zero-inflation

Wrong model

m <- glmmTMB(observedResponse ~ Environment1 + (1|group),
        family = poisson(), data = overData)
res <- simulateResiduals(m)
plot(res)
testZeroInflation(res)

Zero-inflation = 5.16, p-value = 0.

Solution

m <- glmmTMB(observedResponse ~ Environment1 + (1|group),
            ziformula = ~ 1,  # zero-inflation formula
        family = poisson(), data = data) 
res <- simulateResiduals(m)
plotResiduals(res)
testZeroInflation(res)

Zero-inflation = 1, p-value = 1.